home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / legndr.i < prev    next >
Text File  |  1996-02-12  |  3KB  |  102 lines
  1. /*
  2.    LEGNDR.I
  3.    Compute Legendre polynomials and associated Legendre functions.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /* After Numerical Recipes, by Press, et. al., section 6.6.  */
  11.  
  12. func legndr(l,m, x)
  13. /* DOCUMENT legndr(l,m, x)
  14.      return the associated Legendre function Plm(x).  The X may
  15.      be an array (-1<=x<=1), but L and M (0<=M<=L) must be scalar
  16.      values.  For m=0, these are the Legendre polynomials Pl(x).
  17.      Relation of Plm(x) to Pl(x):
  18.        Plm(x) = (-1)^m (1-x^2)^(m/2) d^m/dx^m(Pl(x))
  19.      Relation of Plm(x) to spherical harmonics Ylm:
  20.        Ylm(theta,phi)= sqrt((2*l+1)(l-m)!/(4*pi*(l+m)!)) *
  21.                            Plm(cos(theta)) * exp(1i*m*phi)
  22.    SEE ALSO: ylm_coef
  23.  */
  24. {
  25.   m+= 0;
  26.   l+= 0;
  27.   if (structof(m)!=long || structof(l)!=long ||
  28.       m<0 || m>l) error, "l, m integers with 0<=m<=l for Plm";
  29.  
  30.   if (m > 0) {
  31.     /* sqrt blows up if x out of range */
  32.     pmm= ((m%2? -1:1) * exp(sum(log(indgen(1:2*m:2))))) *
  33.       (sqrt((1.-x)*(1.+x))^m);
  34.   } else {
  35.     if (anyof(abs(x)>1.0)) error, "-1<=x<=1 for Plm";
  36.     pmm= array(1.0, dimsof(x));
  37.   }
  38.  
  39.   if (l==m) return pmm;
  40.  
  41.   pmmp1= (2*m+1)*x*pmm;
  42.   if (l==m+1) return pmmp1;
  43.  
  44.   for (ll=m+2 ; ; ++ll) {
  45.     c1= double(ll-m);
  46.     c= (ll+m-1)/c1;
  47.     c1= (2*ll-1)/c1;
  48.     pll= c1*x*pmmp1 - c*pmm;
  49.     if (ll==l) return pll;
  50.     pmm= pmmp1;
  51.     pmmp1= pll;
  52.   }
  53. }
  54.  
  55. func ylm_coef(l,m)
  56. /* DOCUMENT ylm_coef(l,m)
  57.      return sqrt((2*l+1)(l-m)!/(4*pi*(l+m)!)), the normalization
  58.      coefficient for spherical harmonic Ylm with respect to the
  59.      associated Legendre function Plm.  In this implementation,
  60.      0<=m<=l; use symmetry for m<0, or use sines and cosines
  61.      instead of complex exponentials.  Unlike Plm, array L and M
  62.      arguments are permissible here.
  63.      WARNING: These get combinitorially small with large L and M;
  64.      probably Plm is simultaneously blowing up and should be
  65.      normalized directly in legndr if what you want is Ylm.  But
  66.      I don't feel like working all that out -- if you need large
  67.      L and M results, you should probably be working with some
  68.      sort of asymptotic form anyway...
  69.    SEE ALSO: legndr
  70.  */
  71. {
  72.   m+= 0;
  73.   l+= 0;
  74.   if (structof(m)!=long || structof(l)!=long ||
  75.       anyof(m<0) || anyof(m>l))
  76.     error, "l, m integers with 0<=m<=l; use symmetry for m<0";
  77.  
  78.   lm= l-m;
  79.   not_scalar= dimsof(lm)(1);
  80.   if (not_scalar) result= array(1.0, dimsof(lm));
  81.   else result= [1.0];
  82.   iactive= indgen(numberof(lm));
  83.   for (list=where(lm) ; numberof(list) ; list=where(lm)) {
  84.     lm= lm(list);
  85.     iactive= iactive(list);
  86.     result(iactive)*= lm;
  87.     lm-= 1;
  88.   }
  89.  
  90.   lm= l+m;
  91.   iactive= indgen(numberof(lm));
  92.   for (list=where(lm) ; numberof(list) ; list=where(lm)) {
  93.     lm= lm(list);
  94.     iactive= iactive(list);
  95.     result(iactive)/= lm;
  96.     lm-= 1;
  97.   }
  98.  
  99.   result= sqrt(((2*l+1)/(4.*pi)) * result);
  100.   return not_scalar? result : result(1);
  101. }
  102.